Hydrodynamic crystals: collective dynamics of regular arrays of spherical particles in 

a parallel-wall channel 



M. Baron, J. Blawzdziewicz, and E. Wajnryrjj 

Department 0} Mechanical Engineering, Yale University, P.O. Box 20-8286, New Haven, CT 06520 

(Dated: April 17, 2008) 

Simulations of over 10 3 hydrodynamically coupled solid spheres are performed to investigate 
collective motion of linear trains and regular square arrays of particles suspended in a fluid bounded 
by two parallel walls. Our novel accelerated Stokesian-dynamics algorithm relies on simplifications 
associated with the Hele-Shaw asymptotic far-field form of the flow scattered by the particles. The 
simulations reveal propagation of particle-displacement waves, deformation and rearrangements of 
a particle lattice, propagation of dislocation defects in ordered arrays, and long-lasting coexistence 
of ordered and disordered regions. 



Long-range hydrodynamic interactions between solid 
particles suspended in a fluid result in complex collective 
dynamic phenomena, such as development of ordered ar- 
rays of magnetically driven rotors placed on a liquid in- 
terface [l[ and formation of time-dependent patterns in 
a system of particles immersed in a vibrated fluid Q. 
Collective behavior due to the hydrodynamic coupling 
also occurs in biological systems. A striking example is 
spontaneous formation of vortical arrays of self-propelled 
sperm cells confined to an interface [3(. Hydrodynamic 
coupling also plays an essential role in the synchroniza- 
tion of cilia beating and development of collective waves 
in cilia arrays in small swimming organisms |4|]. 

In confined multiphase systems, the collective particle 
behavior is strongly influenced by bounding walls affect- 
ing the fluid motion. According to recent studies [a, Ifj, |7| , 
hydrodynamic confinement effects are especially signifi- 
cant in parallel- wall channels of width comparable to the 
particle size. Lateral motion of a particle in such a chan- 
nel produces fluid backfiow that is involved in numer- 
ous dynamical phenomena. It enhances relative particle 
motion in confined quasi-2D-suspensions [H, 0], consid- 
erably increases transverse hydrodynamic resistance for 
elongated rigid arrays of spheres moving parallel to the 
channel walls [fjj], and governs propagation of particle- 
displacement waves in linear arrays of drops in a mi- 
crofluidic channel. We show that the fluid backfiow re- 
sulting from particle motion is also responsible for pat- 
tern formation occurring in 2D hydrodynamic crystals 
(i.e. regular particle arrays that are hydrodynamically 
driven) . 

In this Letter we present a numerical study of the dy- 
namics of ID and 2D regular arrays of hydrodynamically 
coupled spherical particles in parallel-wall channels (cf. 
configurations shown in Fig. QJ. We investigate propa- 
gation of displacement waves in linear arrays. In large 
square 2D arrays we report emergence of striking pat- 
terns, such as rearrangements of particle lattice, disloca- 
tion defects, and coexistence of ordered and disordered 
domains. We show that these patterns occur as a result 
of macroscopic deformation of a regular particle lattice, 
and we propose a macroscopic theory describing shape 



evolution of the arrays. 

Our simulations are performed using a novel acceler- 
ated Stokesian-dynamics algorithm to follow evolution of 
about 10 3 particles. Potential applications of our new 
algorithm include studies of collective motion of self- 
propelled particles (e.g. bacterial colonies) in liquid films, 
modeling suspension flows in slit pores, and investiga- 
tions of dynamics of macromolecules (e.g. DNA or poly- 
mer chains) in microfluidic channels. Our acceleration 
technique can also be used in boundary-integral algo- 
rithms for studying dynamics of deformable particles in 
confined geometry. 

Our numerical technique relies on simplifications asso- 
ciated with the far-field asymptotics of the flow scattered 
from the particles. Far from a particle, the scattered flow 
in a parallel-wall channel assumes the Hele-Shaw form, 
i.e. it tends to a 2D parabolic flow that is driven by a 
harmonic pressure distribution Q. In our new approach 
we expand the flow scattered by the particles into a care- 
fully chosen fundamental set of Stokes flows. Close to a 
particle the basis flows form a complete set of solutions 
of Stokes equations in 3D space. In the far-field domain 
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FIG. 1: Particle arrays in parallel- wall channels, (a) System 
definition; (b) lateral displacement wave in linear array; (c) 
square array (top view). 
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FIG. 2: Rescaled dispersion relation for small-amplitude lon- 
gitudinal displacement waves in linear particle arrays in a 
channel of width H = l.ld. Arrays driven by Poiseuille flow 
(dashed lines) and external force (dotted lines) for W/d = 
1.1,1.5,2,3,5,10 (from above). The results for the force- 
driven system are multiplied by the factor a = —0.325. 



p ^> H (where p is the lateral distance from the particle, 
and H is the wall separation) these flows either expo- 
nentially tend to zero or to Hele-Shaw flow driven by a 
2D pressure multipole. The expansion of the flow field 
into the new set of basis fields (obtained by an orthog- 
onal transformation from the fields used in Q) yields 
a sparse system of linear equations, which can be effi- 
ciently solved using iterative sparse-matrix-manipulation 
techniques. Moreover, since the far-field flow is uniquely 
determined by the harmonic pressure distribution, well- 
developed acceleration techniques for Laplace equation 
can be applied to further increase numerical efficiency. 
The simulations discussed below show that our algorithm 
is efficient and highly accurate in both the far-field and 
the near-field regimes. The calculations also reveal sur- 
prisingly rich collective particle dynamics. 

Figures [2] and [3] present our results for propagation of 
particle displacement waves in an infinite train of equally 
spaced particles positioned along a line in the midplane 
of a channel slightly wider than the particle size. The 
particle array is driven either by Poiseuille flow [cf. Fig. 
QJa)] or by a constant external lateral force. We focus on 
the longitudinal waves, where the particle displacements 
8xi from the reference positions Xi = iW (i = 0, 1, 2 . . .) 
on a regular lattice with spacing W occur along the array 
[cf. Fig. HT6)]. 

Figure [5] shows the dispersion relation to = u){k) for 
harmonic displacement waves Sxi = e sin(kxi — ut) in 
arrays with different interparticle spacing. Here k is the 
wave number, u is the wave frequency measured in the 
reference frame moving with the particles, and e < 1 is 
the wave amplitude. The time and frequency are nor- 
malized by the time To in which an isolated particle in a 
channel moves by the diameter d. The shape of the dis- 
persion curves is reflected in the evolution of wave pack- 
ets depicted in Fig. [3l For small interparticle spacing 
the maximum frequency is shifted towards smaller wave 
vectors, because the lubrication forces hinder the relative 
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FIG. 3: Evolution of a wave packet in an array with particle 
spacing W/d = 1.1 (a-c) and W/d = 3 (d-f). Particle dis- 
placements Sxi (i — 0,1, . . .) are normalized by the magnitude 
of the initial perturbation. 



particle motion. Hence, there is a long-wave tail in the 
wave packet shown in Fig. EJ^c). 

In Fig.[2]the frequency u> is plotted rescaled by a factor 
(W/d) 3 to emphasize the universal behavior of the system 
for large values of W. In addition, the results for the 
force-driven train are multiplied by a constant negative 
factor a. We find that for W/d > 5 all rescaled results 
fall onto a single asymptotic master curve. In the regime 
W/d w 2-3 the dispersion relations significantly deviate 
from the master curve, but the scaled results for the flow- 
and force-driven arrays are still nearly identical. Since 
the scale factor a is negative, this behavior indicates that 
for moderate and large interparticle distances the relative 
particle motion in an array driven by an external flow is 
equivalent to the relative motion in a force-driven array 
moving in the opposite direction. 

The above features of the system dynamics result from 
the far-field behavior of the flow field scattered by the 
particles. In the far-field regime an isolated particle sub- 
ject to Poiseuille flow or external force produces the same 
Hele-Shaw flow v HS driven by the two-dimensional dipo- 
lar pressure p HS ~ cos (f>/p (where (j) is the polar an- 
le measured from the direction of the external forcing) 
H|. For W/d > 5 the single-scattering approximation 
corresponding to the superposition of dipolar fields v HS 
adequately describes the system dynamics (and thus the 
rescaled results in Fig. [2] follow a single master curve). In 
the regime W/d w 2-3 evaluation of a single flow reflec- 
tion is insufficient; nevertheless, the results for flow- and 
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FIG. 4: Evolution of a force-driven square array of N — 961 
particles, moving in the diagonal direction. Channel width 
H/d — 1.1 and initial particle spacing W/d — 5 (particles 
are shown magnified by a factor of 2). (a-c) Simulation re- 
sults; (d) prediction of mean- field theory for the same stage 
of evolution as in (a); (e-g) development of a fault line in the 
indicated region. 

force-driven arrays can be rescaled onto each other. This 
is because the first reflection in the multiple-scattering 
sequence for the two systems is nearly identical (apart 
from a rescaling factor), owing to the exponential ap- 
proach of the flow field to the asymptotic Hele-Shaw 
form v HS . With matching initial reflections, the whole 
multiple-scattering sequences for systems with different 
forcing coincide, and the relative particle motion is thus 
the same. This argument is valid not only for linear ar- 
rays but also for other horizontal particle arrangements 
(such as the 2D arrays shown in Figs. [H-[6|). Moreover, 
similar reasoning applies to different kinds of forcing, 
including Marangoni and electrophoretic forces used to 
control particle positions in microfiuidic devices. 

Rich collective phenomena revealed by our simulations 
of 2D hydrodynamic crystals are illustrated in Figs. 0]- 
[5J Figure 2] presents the evolution of a regular square 



array [cf. Fig. QJc)] °f about 10 3 particles undergoing 
diagonal motion produced by a constant force acting on 
all the particles. The initial particle spacing is within the 
far-field asymptotic regime, W/d = 5. Figure [5] shows 
corresponding results for lateral motion of the array. 

Our simulations demonstrate that at short times a de- 
forming square array retains its initial particle ordering 
[cf. Fig. Ufa)]. Subsequently, the system develops some 
striking structural features. Several rows of particle pairs 
separate from the main body of the array, forming a 
shape similar to airplane wings. The front part of the 
array has an approximately hexagonal particle ordering, 
and the middle part retains the square ordering. The 
rear part [marked region in Fig. 0Jc)] has a square par- 
ticle arrangement but with a different orientation than 
the original one. The blowup in Fig. BJe-g) shows that 
the particle rearrangement involves discontinuous parti- 
cle displacements along a "fault line" at the symmetry 
axis of the array. A similar dislocation event (but with- 
out lattice reorientation) is observed in an array in the 
lateral motion (cf. Fig. [5]). There also occur instabili- 
ties responsible for emergence of disordered domains of 
chaotic particle motion. In the diagonal motion [Fig. 
|4j& ; c)] the instability starts at the junction between the 
wings and the body of the array. For the lateral motion, 
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FIG. 5: Same as Fig. [4] but for lateral motion of the array. 
(a, b) Simulation results; (c) prediction of mean-field theory 
for the same stage of evolution as in (a); (d-f) development 
of a dislocation line in the indicated region. 



4 



FIG. 6: Square arrays of (a) N = 961 and (b) N = 256 
particles in diagonal motion. Channel width H/d = 1.1 and 
initial particle spacing W/d = 2. 

the order-disorder transition occurs when the dislocation 
lines approach the center of the array [Fig. HJfe)]. 

Figures [4] and [5] demonstrate that the ordered regions 
can withstand large macroscopic deformations and ran- 
dom perturbations originating from the disordered do- 
mains. Moreover, the ordered crystalline domains can 
rearrange and heal themselves along fault or dislocation 
lines. The strong propensity to maintain the ordered 
structure results from the dipolar hydrodynamic interac- 
tions of neighboring particles. An array also undergoes 
a macroscopic deformation resulting from the combined 
long-range effect of the dipolar flow fields produced by 
individual particles. 

For low-density arrays (i.e. for W/d ^> 1) the macro- 
scopic flow that causes the deformation can be deter- 
mined from the flow field produced by a uniform distri- 
bution of pressure dipoles p HS induced in the array. In 
general, the macroscopic deformation can be described 
using the effective 2D transport equations for suspension 
flow in a parallel-wall channel, u = — t/ p Vp + Vff, and 
j p = —fip'Vp + fXff. Here u is the suspension velocity av- 
eraged across the channel, j p is the particle flux, f denotes 
the density of the lateral force acting on the particles, p 
is the macroscopic pressure, and v a , p a (a — p, f) are 
linear transport coefficients. The suspension velocity and 
suspension flux satisfy the continuity equations V • u = 
and dp/dt = — V • j p , where p is the suspension density. 

The macroscopic deformation of an array, evaluated 
in the uniform-dipolar-moment approximation, is shown 
in Figs. \M,d) and EJc). The results indicate that our 
macroscopic description reproduces the overall shape of 
the arrays for moderate times (i.e. before the complex 
structural features develop). 

We note that the macroscopic equations predict fin- 
gering instabilities near the array corners. In low-density 
arrays (cf. Figs. 2] and [5]) such instabilities are suppressed 
due to the array "stiffness" associated with its ordered 
structure. However, for denser arrays (cf. Fig. [6} the 
macroscopic deforming forces are sufficiently strong to 
destabilize the tips of the array wings, in agreement with 
our macroscopic theory. The results in Fig. [5] indicate 
that the lengthscale for the fingering instability in par- 



ticle arrays is determined by the particle lattice. In our 
effective medium theory there is no intrinsic lengthscale, 
so the size of the fingers in Figs. \M,d) and[5jc) is set by 
the initial condition (i.e. a square with rounded corners). 

The patterns we observe in 2D hydrodynamic crys- 
tals have analogies in other athermal systems. For ex- 
ample, dislocations and chaotic dynamics develop in ar- 
rays of convective cells in a fluid undergoing Benard 
convection [8j and in vibrated granular media [9j. Our 
system has a number of interesting distinctive features. 
First, the pattern formation occurs in the linear Stokes- 
flow regime, and the nonlincarity stems entirely from the 
position-dependence of the multiparticle mobility matrix. 
Next, the dipolar hydrodynamic interactions that main- 
tain particle ordering are non-isotropic (causing, e.g., lat- 
tice reorientation) . Finally, the dipolar flow v not only 
maintains the ordered structure on the local level but also 
produces the macroscopic deformation of the array, lead- 
ing to lattice instabilities. 

Regular particle arrays can be assembled using holo- 
graphic optical tweezers so the collective dynamic 
phenomena revealed by our study should be accessible 
experimentally. The effect of hydrodynamic coupling on 
the motion of regular particle arrays could also be ob- 
served in flow-driven 2D colloidal crystals. The equiv- 
alence of the relative particle motion in systems with 
different forcing can be used to separately control the 
relative particle positions and the position of the center 
of mass of an array. 
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